Initialisation et génération des données
## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
rho2 = 0.1,
mu = c(5,90,5),
omega2 = c(0.5,0.1,0.01),
#Survival data,
nu2 = 0.5,
a = 90,
b = 50,
alpha = 7,
beta = 10)
F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(60,120, length.out = 10) #value of times
# --- longitudinal data --- #
G = 40 ; ng = 8 #nombre de groupe et d'individu par groupe
n <- G*ng*length(t) ; N <- G*ng #nombre total de data, et #nombre total d'individu
dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs
dt_SF <- SF_obs(dt_NLME, param, m)


Statistique exhaustive
\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]
sigma2_a <- sigma2_b <- sigma2_alpha <- 0.1
#Petite fonction pour retourner rapidement l'appel Phi[attr(Phi, i)] où i est '1', '2', ...,
`%a%` <- function(x,var){
if(length(var)== 1) return(x[ attr(x,as.character(var)) ])
lapply(var, function(v) x%a%v) %>% unlist
}
Phi <- fct_vector(function(sigma2, rho2, mu, omega2, nu2, a, b, alpha, beta) {
c(- n/(2*sigma2), -N/(2*rho2), #1, 2
- G/(2*omega2), G*mu/omega2, #3, 4
- G/(2*nu2), #5
a/sigma2_a, -1/(2*sigma2_a), #6, 7
b/sigma2_b , -1/(2*sigma2_b), #8, 9
alpha/sigma2_alpha, -1/(2*sigma2_alpha), #10, 11
beta ) },#12
dim = c(1,1,F., F., rep(1,8)) )
zeta <- function(beta, eta, phi, gamma, a,b, alpha)
{
lbd <- function(t,g) t^{b-1} * exp(alpha*m(t, eta[g], phi[g,]))
P1 <- 1:nrow(dt_SF) %>% sapply(function(i) integrate(function(t) lbd(t,dt_SF$gen[i]) , 0, dt_SF$obs[i])$value )
P2 <- 1:nrow(dt_SF) %>% sapply(function(i) exp(beta*dt_SF$U[i] + gamma[dt_SF$gen[i]] ) )
beta*sum(dt_SF$U) - b*a^-b * sum(P2*P1)
}
zeta.der <- function(beta, eta, phi, gamma, a, b, alpha)
{
lbd <- function(t,g) t^{b-1} * exp(alpha*m(t, eta[g], matrix(phi[g,], nrow = 1)))
P1 <- 1:nrow(dt_SF) %>% sapply(function(i) integrate(function(t) lbd(t,dt_SF$gen[i]) , 0, dt_SF$obs[i])$value )
P2 <- 1:nrow(dt_SF) %>% sapply(function(i) dt_SF$U[i]*exp(beta*dt_SF$U[i] + gamma[dt_SF$gen[i]] ) )
sum(dt_SF$U) - b*a^-b * sum(P2*P1)
}
S <- fct_vector(function(eta, phi, gamma, a, b, alpha) mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ),
function(eta, phi, gamma, a, b, alpha) mean(eta^2), #2
function(eta, phi, gamma, a, b, alpha) apply(phi^2, 2, mean), #3
function(eta, phi, gamma, a, b, alpha) apply(phi, 2, mean), #4
function(eta, phi, gamma, a, b, alpha) mean(gamma^2), #5
function(eta, phi, gamma, a, b, alpha) a, #6
function(eta, phi, gamma, a, b, alpha) a^2, #7
function(eta, phi, gamma, a, b, alpha) b, #8
function(eta, phi, gamma, a, b, alpha) b^2, #9
function(eta, phi, gamma, a, b, alpha) alpha, #10
function(eta, phi, gamma, a, b, alpha) alpha^2, #11
function(eta, phi, gamma, a, b, alpha) #12
uniroot(function(beta)zeta.der(beta, eta, phi, gamma, a,b,alpha), lower = -100, upper = 100)$root,
dim = c(1,1,F., F., rep(1,8)) )
SAEM
initialisation
M <- 1 #nombre de simulation
# --- Initialisation des paramètres --- #
para <- param %>% sapply(function(x) x + runif(1))
para$rho2 = 0.2 ; para$mu <- c(6,3,1) ; para$omega2 <- rep(.1,3)
# --- Initialisation des chaines MC : Z_0 --- #list(attr(dt_NLME,'eta')),#list(attr(dt_NLME,'phi')),#
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, para$rho2) %>% matrix(ncol = 1) ),
phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ),
gamma = list(attr(dt_SF,'gamma')),#1:M %>% lapply(function(i) matrix(rnorm(G, 0, para$nu2), ncol = 1) ),
a = list(matrix(para$a)),
b = list(matrix(para$b)),
alpha = list(matrix(para$alpha)) )
Etape simulation et maximisation du SEAM
sim <- function(niter, Phih, eta, phi, gamma, a, b, alpha)
{
M <- length(phi)
eta <- 1:M %>% lapply( function(i)
MH_High_Dim_para_future(niter, eta[[i]], sd = .036, loglik.eta, phi[[i]], Phih, cores = ncores ))
phi <- 1:M %>% lapply( function(i)
MH_High_Dim_para_future(niter, phi[[i]], sd = c(.028, .036, .021), loglik.phi, eta[[i]], Phih, cores = ncores ))
# gamma <- 1:M %>% lapply( function(i)
# MH_High_Dim_para_future(niter, gamma[[i]], sd = .03, loglik.gamma,
# eta[[i]], phi[[i]], a[[i]], b[[i]], alpha[[i]],
# Phih, cores = 1 ))
a <- 1:M %>% lapply( function(i) MH_High_Dim_para_future(niter, a[[i]], sd = .02, loglik.a,b[[i]], Phih, cores = 1 ))
b <- 1:M %>% lapply( function(i) MH_High_Dim_para_future(niter, b[[i]], sd = .02, loglik.b, Phih, cores = 1 ))
alpha <- 1:M %>% lapply( function(i) MH_High_Dim_para_future(niter, alpha[[i]], sd = .02, loglik.alpha, Phih, cores = 1 ))
list(eta = eta, phi = phi, gamma = gamma, a = a, b = b, alpha = alpha)
}
maxi <- function(S)
{
list(sigma2 = S%a%1,
rho2 = S%a%2,
mu = S%a%4,
omega2 = S%a%3 - (S%a%4)^2,
nu2 = S%a%5,
a = S%a%6, b = S%a%8,
alpha = S%a%11, beta = S%a%12 )
}